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SUPERSONIC PROPULSION SIMULATION BY INCORPORATING COMPONENT 


MODELS IN THE LARGE PERTURBATION INLET 

(LAPIN) COMPUTER CODE 

Gary L. Cole and Jacques C. Richard 
National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135 


SUMMARY 

An approach to simulating the internal flows of supersonic propulsion systems is pre- 
sented. The approach is based on a fairly simple modification of the publicly-available Large 
Perturbation Inlet (LAPIN) computer code. LAPIN uses a quasi-one-dimensional, inviscid, 
unsteady formulation of the continuity, momentum, and energy equations. The equations are 
solved using a shock-capturing, finite-difference algorithm. The original code, developed for 
simulating supersonic inlets, includes engineering models of unstart/restart, bleed, bypass, and 
variable duct geometry, by means of source terms in the equations. The source terms also pro- 
vide a mechanism for incorporating, with the inlet, propulsion system components such as com- 
pressor stages, combustors, and turbine stages. This requires each component to be distributed 
axially over a number of grid points. Because of the distributed nature of such components, this 
representation should be more accurate than a lumped-parameter model. Components can be 
modeled by performance map(s), which in turn are used to compute the source terms. The 
report describes the general approach first. Then, simulation of a compressor /fan stage is 
discussed to demonstrate the approach in detail. Generalized source terms, valid for either 
subsonic or supersonic flow through the compressor, are derived for the momentum and energy 
equations. A starting transient is presented as a demonstration case for a single-stage, super- 
sonic fan in a mixed-compression inlet. Discussion of the compressor stage implementation in 
LAPIN provides guidance for the addition of other components. 


INTRODUCTION 

A computer simulation can be a valuable tool during a propulsion system development 
cycle. In particular, dynamic simulations provide a means for investigating potential inter- 
actions between individual components. For example, supersonic propulsion systems may 
exhibit an inlet unstart (expulsion of the normal shock) that can induce engine stall (refs. 1 to 3) 
and possibly surge. Conversely, engine stall will cause inlet unstart (ref. 2). Simulations can 
also be used to guide the design and testing of propulsion control systems prior to wind tunnel 
and flight tests (refs. 4 and 5). 

Typically, lumped-parameter models (e.g., refs. 6 and 7) have been used to represent tur- 
bine type engines. More recently, a distributed turbojet engine model, using a one-dimensional 
finite-difference approach, was proposed in reference 8. Distributed models are generally more 
accurate, while lumped models have the advantage of faster execution times. Lumped models 
(e.g., ref. 9) have also been used to represent the gas dynamics of supersonic inlets. However, 
the ability to represent normal shocks with such models is quite limited. A more accurate 
model, using a linearized wave-equation approach, is described in references 10 and 11. That 



model is able to predict inlet resonance conditions with reasonable accuracy, but is limited to 
small perturbations or excursions in normal-shock position. 

The Large Perturbation Inlet (LAPIN) computer code was developed (ref. 12) to deal 
with large flow-field perturbations in supersonic inlets, such as hammer shock and unstart/ 
restart transients. The approach is based upon a quasi-one-dimensional, inviscid, unsteady 
formulation of the continuity, momentum, and energy equations. The equations are solved using 
one of several user-selectable shock-capturing, finite-difference algorithms. LAPIN is publicly 
available from the NASA Computer Software Management and Information Center (COSMIC). 

In many cases, propulsion system simulation can be adequately accomplished by coupling 
a distributed inlet model with a lumped engine model. One exception is a Supersonic Through- 
Flow Fan (STFF) system. The STFF system exhibits a so-called "starting transient” — the pas- 
sage of a normal shock through the supersonic fan stage (from upstream to downstream) during 
the acceleration phase of the flight path. Unstarting, the reverse transient, occurs during 
deceleration. A STFF system configuration proposed for commercial supersonic transports and 
its advantages relative to a conventional mixed-flow turbofan engine are described in refer- 
ence 13. The simulation of a STFF system and its starting/ unstarting transients requires the 
use of distributed inlet/engine models like LAPIN and the engine model reported in reference 8. 
The simulation effort described in this report was motivated by NASA Lewis’ interest in the 
STFF concept. However, the approach is applicable to conventional supersonic propulsion sys- 
tems as well. 

The simulation approach, described in this report, is based on a modified version of the 
LAPIN code. LAPIN was originally developed for the simulation of mixed-compression, super- 
sonic inlets. The effects of an engine at the diffuser exit had to be approximated by appropriate 
variations of the downstream boundary conditions. Engineering models of unstart/restart, 
bleed, bypass, and variable duct geometry effects were implemented by means of source terms in 
the governing equations. Subsequent studies have shown that the entire compression system 
(inlet and fan/compressor stages) can be simulated by adding source terms in the LAPIN model 
and by using component performance maps to furnish the information needed to calculate the 
additional source terms. This approach is generic and can be extended to incorporate other 
propulsion components with the inlet, including combustors and turbine stages. The turbine- 
engine model described in reference 8 uses a similar approach. 

This report describes the general approach for adding propulsion-component source terms 
to LAPIN and the terms specific to fan (compression) stages. A starting transient is presented 
as a demonstration case for a single-stage, supersonic fan in a mixed-compression inlet. The 
actual implementation and modifications to LAPIN subroutines is discussed, which provides 
guidance for the addition of other components. 


GENERAL APPROACH 

The general approach for adding propulsion-component source terms to the LAPIN equa- 
tions is illustrated in figure 1. The entire flow field is represented by a one-dimensional grid 
with supersonic boundary conditions (BCs) upstream and subsonic boundary conditions down- 
stream. LAPIN solves the unsteady, inviscid continuity, momentum, and energy equations 
(fig. 1) in finite-difference form at each grid point. The equations are written in weak conserva- 
tion law form, with source/sink terms on the right-hand side. The pressure times duct-area 
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axial gradient and pressure times duct-area time-rate-of-change terms are included with the 
momentum and energy source terms, respectively. The quantities M g , F g , and Q g represent 
the remaining cumulative source terms for the continuity, momentum, and energy equations, 
respectively. They are used to incorporate: terms originally included for mass bleed/bypass and 
unstart/restart phenomena; terms added later (ref. 14) for mass injection, duct pressure recovery 
losses, and heat addition; and the new terms proposed in this report for incorporating propulsion 
cornponents. 

In general, a component can be located anywhere on the grid and can be modeled by some 
form of performance map. (A component model in equation form could also be used.) The 
source terms must be derived in terms of information that can be furnished by the map(s) or 
equation(s). In LAPIN, the source terms are calculated at each time-step using the following 
procedure. First, inputs to the map are derived from a combination of flow field information 
and "externally” supplied information. Examples of external information are mechanical speed 
for a compressor/turbine component and fuel flow for a combustor. The external information is 
supplied either by means of LAPIN input statements (NAMELISTs) or an auxiliary subroutine 
(e.g., a control system algorithm). Outputs from the performance map, needed for the source 
terms, are then computed. Finally, component source terms required for the continuity, energy, 
and momentum equations are calculated and distributed to the grid points encompassed by the 
component. 

To simulate a particular "engine,” the associated components (e.g., compressor, com- 
bustor, turbine) are located in tandem on the grid. It is conceivable that component source 
terms will overlap other source terms. An example might be a compressor with a bleed region. 
The next section discusses some of the details of adding source terms to LAPIN that are specific 
to compression-stage components. 


COMPRESSION STAGE SOURCE TERMS 

The performance of a compressor or fan component is typically represented by a multi- 
dimensional map like the one illustrated in figure 2. In the illustration, component total 
pressure ratio and adiabatic efficiency are plotted as functions of corrected mass flow rate and 
corrected rotor speed. Such a map typically represents the overall performance of the compo- 
nent. However, particular applications and/or desired levels of simulation fidelity may require 
separate maps for each stage (i.e., rotor-stator combination) or even separate maps for each 
rotor and stator. Another approach might represent the rotor by a map and the stator by a 
simple drag source term. For the map shown in figure 2, the only external input required is 
compressor (or fan) mechanical speed. Other variables needed to compute the map inputs are 
mass flow rate and the total temperature and total pressure at the stage inlet. The map outputs 
are the stage total pressure ratio and adiabatic efficiency. As an alternative, the map inputs 
could be defined as corrected speed and stage total pressure ratio. In this case, the map outputs 
would be the stage inlet corrected mass flow rate and the adiabatic efficiency. In some cases, the 
shape of the corrected speed lines will dictate which procedure is used. 

For component maps, such as the one shown in figure 2, it is possible to have certain com- 
binations of inputs (i.e., corrected mass flow rate and corrected speed) that extrapolate to a map 
output that lies outside the normal operating range defined by the map. A computational strat- 
egy must be devised to handle such cases. If an output point should occur in the stall region, 
two possible options would be to (1) output a message that stall has occurred and terminate the 
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simulation or (2) use an estimated stall behavior (i.e., reduced total pressure ratio and adiabatic 
efficiency) and continue the simulation. In most cases, the second choice is preferable. That 
approach was adopted in reference 8 for stage characteristics like those in figure 3. The curves 
drawn to the left of point b in figure 3 are assumed to represent the in-stall performance. Refer- 
ences 15 and 16 present the so-called <f> versus ^ characteristics, which are similar to the 
pressure ratio versus corrected mass flow rate curve of figure 3(a). The </> and ^ parameters 
are normalized with respect to the mean rotor speed. Both ^ an d corrected mass flow rate are 
indicators of stage inlet axial velocity or Mach number. Reference 15 explains that the <j> 
versus if) curves remain the same, with varying rotor speed, if the compressibility effects are 
small. Hence, the single-curve characteristic is valid at low Mach numbers and may provide a 
useful first approximation at moderately high Mach numbers. This would not apply to a super- 
sonic stage and it is obvious that the approach and type of performance map(s) selected to simu- 
late compression stages will be highly dependent on the particular application. Additionally, a 
time lag may be imposed on map inputs/outputs, as suggested in reference 8. 

For purposes of illustrating the source terms that are specific to compression stages, it is 
assumed that the inputs to the performance map are corrected speed and stage inlet corrected 
mass flow rate. The corresponding performance map outputs are the stage total pressure ratio 
and adiabatic efficiency, which can be combined to calculate the stage total temperature ratio. 

Since compression stages add energy and momentum to the flow, the corresponding source 
terms must be derived for the energy and momentum equations. Details of the derivation are 
presented in appendix B. The sketch in appendix B shows a representative compression stage 
and the associated nomenclature. The equations for the source terms are as follows: 


Energy Equation 

The total work per pound done by the stage is equal to the change in enthalpy across the 
stage (ref. 17). 


Ah - c p ^T t3 - T t2 ) 

Knowing the stage inlet total temperature, pressure, and mass flow rate and the stage mechani- 
cal speed, the stage total pressure ratio and adiabatic efficiency can be obtained from a perform- 
ance map (e.g., fig. 2). The change in total temperature across the stage (ref. 17) is then 
calculated from: 


T t3 " T t2 ~ T t2 


Pt3 

(i-i )h 

- l 

l P t2, 



The stage total rate of energy addition, written in terms of known flow field information and 
information derived from the component performance map, is 
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Qsfc - m 2 c p( F t3 " T t2) “ (M u )2 c p(Tt3 ” T fc2 ) 


Momentum Equation 

Between two stations, i and i+1, within the compression stage region, the net force due 
to the fan in the axial direction can be expressed as: 

F st = m i+l u i+l - m i u i + A i+lPi+l - A iPi - Pm( A i+l - A i) 

This is a general expression that accounts for flow-area variation through the compression stage 
region. As shown in appendix B, the force can be rewritten as 


F„ - ftufA.FAC, ♦ p,A l FAc J A|+ '" A ‘) 

2 A; 


which is in terms of known flow field variables at the cell entrance (station I) and the informa- 
tion derived from the component performance map (total temperature and pressure ratios across 
the cell contained in FAC 1 and FAC 2 ). This term represents only the incremental force 
between stations i and i+1. The total force term for the entire stage cannot be represented by 
a single equation, when flow area varies. If the stage flow area is constant, then the previous 
equation reduces to: 


F st - P 2 n 2 A 2 FAC x + p 2 A 2 FAC 2 

which is the total force between the stage endpoints (stations 2 and 3). FACj and FAC 2 are 
the same terms from the previous equation, evaluated at the stage endpoints (i.e., i+1 is the 
stage exit 3, and i is the stage inlet 2). 

The governing equations (fig. 1), as solved by LAPIN at each grid point, are normalized 
and in finite difference form. In general, the source term contribution at a grid point covers a 
cell region that extends between stations (i) and (i+1) that are midway between the grid point 
(j) and the adjacent upstream (j-1) and downstream (j+1) grid points. Cell lengths depend on 
the grid spacing and grid points at the component end points are handled specially. The nor- 
malized source terms per unit length, specific to compression stages, are derived in appendix B 
and shown to be as follows: 
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Energy Equation 


Qst = W 2 


; t 43 - T t2 i 

T ^ 

L ref 

k Tt ° , 

jstage , 


which is the total energy source term to be distributed to the compression stage grid points. 


Momentum Equation (Variable Area Through Stage) 


^st “ 


(mf 


pi 


FAGj + pjFAC 2 


( A i+1 + A i) 


2 A: 


J ref 


x i+i - x i 


which is the incremental momentum source term between stations i and i+1. 


Momentum Equation (Constant Area Through Stage) 


** = 


(Wif 


P 2 


FAC X + p 2 FAC 2 


J ref 


J sfcage 


which is the total momentum source term to be distributed to the compression stage grid points. 
FAC X and FAC 2 are the same terms from the previous equation, evaluated at the stage end- 
: points (i.e., i+1 is the stage exit 3, and i is the stage inlet 2), and L atage = x 3 - x 2 . 


DEMONSTRATION CASE 

In order to demonstrate the actual implementation of compression stage source terms, an 
example case was contrived and run using LAPIN. The example case corresponds to the start- 
ing transient for a single-stage supersonic fan in a mixed-compression inlet. The compression 
stage was implemented by modifying the heat addition section of the LAPIN subroutine USS. 
The implementation of that source term is explained in reference 14. Details of the required 
FORTRAN modifications, including the addition of a new NAMELIST and changes to sub- 
routines INPUT, BBINIT, and USS are addressed in appendix C. 


Description 

This section describes how the demonstration transient was simulated and the required 
LAPIN input, A listing of the input and results for selected time steps is given in figure 4. 
LAPIN has seven NAMELIST inputs, each beginning with &NAME and ending with &END. 
NAMELIST CONTRL is used to specify the finite-difference algorithm (split characteristics for 
this case, ref. 12), the time step size, and the time steps to output results. Boundary condition 
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initial values and time histories are defined by NAMELIST BC. The free-stream Mach number, 
static pressure, and static temperature initial conditions are MHUSIC (MacH UpStream Initial 
Condition), PRUSIC, and TMUSIC respectively. The initial exit Mach number MHDSIC is 
0.164 and is varied with time, as will be explained later. NAMELIST GEOM is used to define 
the duct geometry. A geometry, without a centerbody, was chosen for simplicity. In an actual 
case a centerbody is required for the fan and would be used to control duct Mach number dis- 
tribution. The duct axial-area distribution is shown in figure 5. Distances are normalized by 
the cowl radius, COWLR. The duct throat is located at x = 0.7 and has a Mach number of 
about 1.5 for the specified free-stream conditions. As shown in figure 5, the fan location (spec- 
ified in NAMELIST HEAT) has been set in a region where the duct area varies, to demonstrate 
that the source terms are not restricted to constant area sections. The duct was configured to 
give a fan face Mach number of about 2.0, when the fan is operating supersonically. Although 
not used in the demonstration case, fan-stage blockage could be specified in NAMELIST GEOM 
by using the strut area and location arrays, ASTRUT and XSTRUT respectively. Bypass, bleed 
and injection mass flows are not used in the demonstration case, but could be defined in 
NAMELIST BYBL and NAMELIST INJECT if required. NAMELIST HEAT is ordinarily 
used to specify location and rate of heat addition to the flow. As mentioned earlier, the heat 
addition section of LAPIN subroutine USS was modified to simulate the effects of a compression 
stage. Values of stage adiabatic efficiency and total pressure ratio were input by means of 
NAMELIST HEAT, rather than using a map subroutine to determine component performance. 
A single compression stage (NHEAT = 1) is used, with normalized starting location XHEAT(l) 
of 1.56 and length LHEAT(l) of 0.7 (fig. 5). The stage adiabatic efficiency was assumed to 
have a constant value of 0.65 for all values of pressure ratio and was specified by using 
XHEAT(2). The HEATIN array defines linear variations of stage total pressure ratio between 
the time points specified by the TIMHT array. NAMELIST LOSS is used to define drag and 
skin friction terms, not used by the demonstration case. 


Discussion of Results 

The duct exit Mach number and stage total pressure ratio were varied in a manner 
intended to simulate a supersonic fan starting transient. Their time histories are plotted in 
figure 6. Axial distributions of Mach number, static pressure, total pressure, and total tempera- 
ture are shown in figure 7 at several instants of time during the transient. Figure 8 shows time 
histories of the stage total pressure ratio and temperature ratio during the transient. The ratios, 
initially 1.0, were calculated from downstream (x = 2.4) and upstream (x = 1.5) values. 

During the first 0.03 sec of the transient, the simulation is brought to a steady-state con- 
dition. Curve A of figure 7 shows the initial operating condition with the normal shock located 
upstream of the fan stage and just downstream of the inlet throat (x = 0.7). The normal shock 
location is evidenced by the sharp decrease in Mach number (fig. 7(a)) and increase in static 
pressure (fig. 7(b)). The normal shock also causes a drop in total pressure that is evident in 
curve A of figure 7(c) for x greater than about 0.8. It should be noted that the same curve 
also exhibits a total pressure increase in the vicinity of the normal shock. This increase is the 
result of a numerical instability that is often produced by finite-difference algorithms at grid 
points adjacent to discontinuities like normal shocks. The severity of the instability depends on 
the grid spacing and the Mach number upstream of the shock. A more severe case is observed in 
curve E of figure 7(d) at about x = 3.1. 
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With the fan initially operating at a subsonic Mach number, the fan total pressure ratio is 
ramped from 1.0 to 1.4 during the timespan from 0.03 to 0.07 sec (fig. 6(b)). The normal shock 
moves downstream, to satisfy the duct exit Mach condition, as the fan pressure ratio rises. As a 
result, the fan inlet Mach number varies. The shock position, for the pressure ratio of 1.4, is 
just upstream of the fan (fig. 7, curve B). The expected total temperature ratio is 1.155 based 
on a total pressure ratio of 1.4 and an efficiency of 0.65. Figure 8 verifies that the LAPIN 
energy and momentum source terms do produce the correct stage total pressure and temperature 
ratios. 

From 0.09 to 0.12 sec, the fan pressure ratio remains constant while the duct exit Mach 
number is ramped up to "pop" the shock through the fan. Once the shock enters the fan, the 
fan inlet Mach number jumps to a supersonic Mach number of about 2.0 and remains constant. 
At 0.1 sec, the shock is approximately halfway through the fan (fig. 7, curve C). The fan overall 
total pressure ratio momentarily drops below 1.4 with the shock in the fan, due to total pressure 
losses associated with the normal shock. As the shock moves downstream, the shock Mach num- 
ber increases and the associated shock total-pressure loss increases. The resulting decrease of the 
fan overall pressure ratio is illustrated in figure 8(a). The spikes in pressure ratio, at values 
above 1.4 and below 0.85 in figure 8(a), are due to the aforementioned numerical instability. 
After the shock exits downstream of the fan, the pressure ratio jumps back to a value of 1.4 
(fig. 8(a)). Curve D of figure 7 shows conditions at t = 0.11 sec, after the shock exits the fan 
but before completion of the ramp in duct exit Mach number. 

The final phase of the transient results from an increase in fan pressure ratio from 1.4 to 
an assumed operating value of 2.4. The final steady-state conditions are shown by curve E of 
figure 7. Again, figure 8 shows that the LAPIN energy and momentum source terms do produce 
the correct stage total pressure and temperature ratios. The expected total temperature ratio, 
corresponding to a total pressure ratio of 2.4 and an efficiency of 0.65, is 1.437. 

The authors recognize that the demonstration transient does depart from a realistic 
starting transient in some respects. However, the main objective of the transient was to 
demonstrate that the source-term approach for representing components does work. In par- 
ticular, the compression stage source terms were shown to be valid for both subsonic and 
supersonic conditions at the stage inlet. 

CONCLUDING REMARKS 

An approach for simulating supersonic propulsion systems, using the publicly-available 
Large Perturbation Inlet (LAPIN) computer code, was described. In general, source terms in 
the governing continuity, energy, and momentum equations can be used to incorporate compo- 
nent models with a supersonic inlet. Simulation of a compression/fan stage was discussed to 
demonstrate the approach in detail. The source terms, required for the momentum and energy 
equations, were derived in terms of quantities that can be obtained from a compression stage 
performance map. The source terms are generic and are valid for stages with either subsonic or 
supersonic inlet conditions and with variable flow area through the stage. Results were pre- 
sented to demonstrate the approach for the starting transient of a supersonic fan in a mixed- 
compression inlet. The results verified that the energy and momentum source terms do produce 
the correct stage total pressure and temperature ratios. Details concerning the LAPIN imple- 
mentation of the source terms were presented, providing guidance for incorporating other com- 
ponents such as combustors and turbines. This will allow LAPIN to be used to simulate a 
variety of supersonic propulsion configurations. 
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APPENDIX A. - SYMBOLS 


A local area 

c p specific heat at constant pressure 

E total internal energy 

F g cumulative momentum source term 

F gt compression-stage momentum source term 

Ah change in enthalpy 

L length 

M Mach number 

M g cumulative continuity source term 

m mass flow rate 

p pressure 

Q g cumulative energy source term 

Q gt compression-stage energy source term 

T temperature 

t time 

u axial fluid velocity 

x axial coordinate 

7 ratio of specific heats 

tj stage adiabatic efficiency 

p fluid density 

Subscripts: 

1 axial section within component region 

j grid point 

m mean value 

ref reference value 

t total or stagnation condition 

0 free-stream condition 

2 compression stage inlet 

3 compression stage exit 
Superscripts: 

normalized variable 
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APPENDIX B - DERIVATION OF COMPRESSION-STAGE 
ENERGY/MOMENTUM SOURCE TERMS 
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Energy Equation 


The total work per pound done by the stage is equal to the change in enthalpy across the 
stage (ref. 17). 

Ah = c p (T t3 - T t2 ) 

The total rate of energy addition of the stage per unit length is 

A u ?2-A-2 u 2 c p \Tt3 “ ^ta) 

Qst = m 2 Ah = 

^stage 


In LAPIN the energy source terms are normalized, as follows 


Qst L re f _ / , 2 A 2 u 2 c p(^t3 " ^2) L ref 

(pU 3 A) ref (pAu) re f u^ ef ^stage 


where u ref is the free-stream sonic velocity. 

The LAPIN primitive state variables are 


so 


pA * u 
— - — - , u = , 

(^A)ref u ref 


and p = 

(pu 2 A) ref 


Qst 


= f>2*2 c p 



L ref 

^stage 


Assuming adiabatic and thermally perfect conditions, free-stream velocity is related to free- 
stream sonic velocity by (ref. 18): 


/ > 
u 0 


r A 

5 

6 

w 

5 

* kJ 

5 
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2 

Rearranging, substituting u 0 /7RT 0 , and multiplying by T t o/ T to 


2 5 + M 0 7RT 0 


U ref = 


6 T, 


L to 


to 


since 


c p - 7 R/(i - 1 ) aid T 0 /T to -5/(5 * M*) 


a _ 5 * M 0 

*ref c 


(7 " l)c c 


5 +M! 


L t0 


which, for 7 = 1 . 4 , reduces to 


U ref =V 


L t0 


2 A 

substituting for u fef , Q st becomes 


A A A 

Qst = ^2 U 2 


3 (^t3 “ , ^'t2 ) ^ref 


T t0 ^stage 


which is the total energy source term to be distributed over the compression stage grid points. 


Momentum Equation 

Between two stations, i and i+1, within the compression stage region, the net force due 
to the fan can be expressed as: 

F st = m i+l u i+l “ m i u i + A i+lPi+l " A iPi " Pm( A i+l - A i) 


This is a general expression that accounts for area variation through the compression region. 
Assuming the mass flow rate is constant (i.e., m i+1 = mj = PjUjA.) and 

Pm = (Pi+i + Pi )/ 2 
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F st = Pi u i A i 


/ > 
Ui+1 - 1 

+ 

++1 H- Aj' 

P;A: 

f \ 

Pi+1 t 

u i 

\ 1 ) 


2Aj 

< 1 J 

"1 '1 

I'p. J 


Using relationships from reference 18 , the ratio Uj +1 /u- can be converted to a function of Mach 
numbers and total temperatures — 


u i+l = u i+l 

7 RTj 


7 RT i+1 


T t 

l i+X 

u i u i -v 

7 RT i+1 

% 

7 RTj % 

Tt i+i % 

T t 


substituting M for U/V7RT and (1 + 0 . 2 M 2 ) 1 / 2 for (T t /T) x / 2 



f \ 

1/2 


u i+l _ Mi+i 

1 + 0 . 2 M 2 
1 

T ‘,„ 

Ui ^ 

1 + 0.2M 2 , 
1+1J 

\ 

T t> 


Likewise, the ratio Pj +1 /pj can be converted to a function of Mach numbers and total 
pressures — 


Pi-hl _ Pi-|-1 P*i+i 
Pi Pi Pt i+1 Pt. 


1 + o.2m: 


2 


Y 5.5 


1 + o.2m: 


2 


i+l 


Pt 


i+l 


Pt, 


Letting MR = (l + 0.2M*) j{\ + 0 . 2 Mj +1 ), F gt becomes 


r 


FACj 




F st = Pju. Aj 


Mi+1 MR 1 / 2 

T ‘i+i j 

+ 

Ai+i + A s 

Pi A i 

MR 3,5 Pti+1 - 1 

[* 4 

T *. J 


2A i 

k 1 J 

Pt > 


FAC, 


The normalized force term per unit length is 


Fst - 


F st L ref 


*1+1 " x 


1 U Z A) ref 
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In terms of LAPIN primitive variables 


Fst = 


feujf 

IfjLFACj + 

r A i+1 + A;' 

Pi FAC, 

L ref 

Pi 

2 A: 

1 J 

rj * * 

x i+l " x i 


which is the incremental momentum source term between stations i and i+1, for the variable 
area case. If area through the fan stage is constant, the expression for F gt becomes 


Fst - 




P 2 


FAC i + P2FAC2 


J ref 


J 8tage 


which is the total momentum source term to be distributed over the compression stage grid 
points and where FACj and FAC 2 are evaluated at the stage end points (i.e., station i+1 is 
the stage exit 3 and station i is the stage inlet 2) and L stage == x 3 - x 2 . 
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APPENDIX C. - LAPIN IMPLEMENTATION OF SOURCE TERMS 

The continuity, momentum, and energy equation source terms provide the hooks for 
incorporating propulsion component models with LAPIN. This appendix begins with a general 
discussion of the procedure used to implement additional source terms, and modifications 
required to some LAPIN computer-program subroutines. Then, implementation of the compres- 
sion stage source terms in subroutine USS is discussed in some detail, in order to provide 
guidance for adding other components. The discussion assumes that the reader has some famil- 
iarity with LAPIN. Additional information concerning the program, its subroutines, and the 
source terms can be found in references 12, 14, and 19. (To the authors’ knowledge, there are no 
references that discuss the LAPIN FORTRAN source code.) 

GENERAL PROCEDURE 

The following steps outline the general procedure for implementing additional source 

terms: 

1. Represent component(s) by means of performance map(s) (e.g., compression stage pressure 
ratio and efficiency as functions of corrected airflow and corrected speed). Three ways to 
include component performance map data in LAPIN are: calls to a subroutine, interpola- 
tion of block data arrays, or interpolation of arrays input by means of NAMELIST(s). 

2. Derive the normalized source terms for the components in terms of quantities that can be 
obtained from the performance maps (e.g., total pressure ratio and temperature ratio for 
compression stages). Source terms for the compression stage example presented in this 
report are derived in appendix B. 

3. Modify LAPIN to reflect the addition of the new source terms to the continuity, momen- 
tum, and energy equations. In general, this requires modification of LAPIN subroutines 
INPUT, BBINIT, and USS, as described below. 

General FORTRAN Modifications and Additions 

Subroutine INPUT . — Subroutine INPUT reads the input data for LAPIN by means of 
NAMELISTs and sets up a number of variables that are used as flags throughout the program. 
New NAMELIST(s) must be added to LAPIN to provide for definition of additional compo- 
nent^) and external input time histories (e.g., mechanical speed as a function of time for 
compression stages). The new NAMELIST(s) can be modeled after the existing ones used for 
other user-supplied source terms — BYBL for bypass and bleed mass flows; INJECT for injection 
mass flow; HEAT for heat addition; and LOSS for drag and skin friction. Subroutine INPUT 
must be modified to include the NAMELIST(s) declaration(s) and the corresponding READ 
statement(s). COMMON block(s) must also be added to pass the necessary information to 
subroutines BBINIT and USS. 

Subroutine BBINIT . — This subroutine is used to determine the starting and ending grid 
indices associated with each user-defined source-term zone. Subroutine BBINIT must be 
modified to include the COMMON blocks added to subroutine INPUT, and statements that 
determine starting and ending grid indices for the new source terms. 
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Subroutine USS . — Subroutine USS is used to calculate the values of the source terms at 
each time step for the continuity, energy, and momentum equations. A separate section must be 
added to subroutine USS to calculate the source terms for each new component type. Compo- 
nents of the same type are calculated in the same section by using arrays of dimensioned vari- 
ables. In general, the calculation procedure is as follows, for each component of the same type: 

1. Set the starting and ending grid indices for the component. 

2. Calculate flow field conditions required as inputs to the component performance map. (For 
the compression stage these are the conditions at the grid point just upstream of the stage, 
as will be explained in more detail below.) 

3. Calculate the external variables required as inputs to the component performance map. 

(For the compression stage this would be a function of the mechanical speed obtained from 
NAMELIST input or from an auxiliary subroutine.) 

4. Calculate the component performance map output variables, needed to calculate the source 
term values. This could be done by calling a subroutine. 

5. Calculate the source term value(s) at each grid point associated with the component. 

6. Add the source term value(s) to the right hand side of the continuity, energy, and momen- 
tum equations. 


MODIFICATIONS TO SUBROUTINE USS FOR COMPRESSION STAGE SOURCE TERMS 

Implementation of the compression stage source terms is illustrated in the FORTRAN list- 
ings of USS provided at the end of this appendix. A listing for the variable flow area case is 
given first, followed by a listing for the constant area case. The heat addition section of USS 
was modified to implement and test the compression stage source terms derived in appendix B. 
For the following discussion, the reader should refer to the first listing. Throughout the dis- 
cussion, FORTRAN program variables are expressed in all capital letters. Corresponding quan- 
tities from appendix B are given in square brackets. 

At line 5, if there are no stages (i.e., NHEAT< = 0) this section of code is bypassed by 
jumping to statement 800. Most of the lines from 6 to 17 are standard for all source term 
sections of USS, JBEG and JEND (lines 7 and 8) are the beginning and ending grid indices for 
the stage. The variable JSTAGE was added to identify the grid point just upstream of the stage 
that will not have a contribution from the stage source term. Flow variables at JSTAGE are 
needed to compute inputs to the compression stage map. As mentioned in the discussion of the 
demonstration transient, a stage map was not used. Instead, the pressure ratio was input as a 
function of time by means of the NAMELIST HEAT arrays HEATIN and TEMHT. The instan- 
taneous stage pressure ratio PR [P t3 /P t2 ] is determined by interpolation using subroutine 
LININT (lines 14 to 19). The adiabatic efficiency ETA [17] is input by means of NAMELIST 
HEAT variable XHEAT(2) on line 20. When a compression stage performance map is used the 
following steps should be taken: 
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1. Calculate corrected airflow at grid point JSTAGE. 

2. Determine mechanical speed (e.g., by interpolation of external input via NAMELIST) and 
compute the associated corrected speed. 

3. Call the performance map subroutine to obtain stage total pressure ratio and adiabatic 
efficiency. 

Once PR and ETA are known, the ratio DTOVRT [(T t3 - T t2 )/T t2 ] is calculated on line 30. 

The variables Q(J,1) [/>], Q(J,2) [pu], and Q(J,3) [E] are the LAPIN primitive variables at 
grid point J (grid point JSTAGE on lines 31 and 32). They are used to calculate TTJSTAGE 
[T t2 /T t0 ] on line 34. The variable CELHRT3 [3(T t3 - T t2 )/T t0 ], needed for the energy source 
term, is calculated on line 36. The variable DPOVRP [(P fc3 - P t2 )/P t0 ], needed for the momen- 
tum source term calculation, is calculated in a similar manner (lines 37 to 42). 

The next step is to distribute the source terms to the grid points encompassed by the stage. 
In general, the source term contribution at a grid point covers a cell region that extends between 
stations (i) and (i+1) that are midway between the grid point (j) and the adjacent upstream 
(j-1) and downstream (j+1) grid points (see sketch in appendix B). Cell lengths depend on the 
grid spacing. Grid points at the component end points are handled specially. 

The loop for calculating the source term contributions at each component grid point begins 
on line 44. The cell length DXCEL [(x }+1 - Xj)/L ref ] for the grid point is calculated first in the 
section from line 45 to 69. This section is common to most of the LAPIN source terms. The 
portion of the total temperature change for the cell is calculated on line 73, based on cell length. 
Finally, the energy source term due to the stage is calculated and added to the cumulative 
energy source term for grid point J, G(J,3), on line 74. PSX(J) is a coordinate transformation 
term. 

Calculation of the momentum equation source term for the grid point begins on line 76. In 
order to calculate this term it is necessary to determine the total pressure ratio and total tem- 
perature ratio and the exit Mach number for each grid cell. The cell inlet total pressure PTXX 
[P t .], total temperature TTXX [T t .], and Mach number XXM [Mj] are based on information from 
the previous time step. These values are obtained from average flow variables, as shown on 
lines 85 to 95. A slightly different procedure is used (lines 99 to 110) for the first grid point 
(JBEG). For this implementation of the compression stage source terms, it is assumed that the 
change in total pressure and temperature across each cell is proportional to the cell length. That 
fact, along with the previously determined values of DTOVRT and DPOVRP, is used to calcu- 
late the cell total temperature and pressure ratios, TTR [T t /T t< ] and PTR {Pj /P t .], on 

lines 115 and 116. The cell exit Mach number XXM2 [M i+1 ] still remains to be calculated. This 
is done using an iterative procedure on lines 117 to 140, where the inlet and exit flow areas for 
the cell were obtained on lines 80 to 83. Note that the iterative procedure depends on whether 
the cell inlet Mach number is subsonic or supersonic. All variables necessary to calculate the 
momentum source term for the compression stage are now available. The term is already pro- 
portioned to the length of the cell by virtue of the method used to calculate the cell total tem- 
perature and pressure ratios. The factors FAC1 [FACj] and FAC2 [FAC 2 ] for the momentum 
source term, are calculated on lines 142 and 143, respectively. Finally, the source term itself is 
calculated on lines 144 to 146 and added to the cumulative momentum source term for grid 
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point J, G(J,2), on line 147. The terms 92 AVG, Q1AVG, and GM1*(Q3AVG - 0.5 
*Q2AVG*Q2AVG/Q1AVG) represent /JjUj , />j , and p i5 respectively. PSX(J) is again a 
coordinate transformation term. 

Calculation of the momentum term at each grid point is simplified, when the flow area 
through the compression stage is constant. The Fortran listing for the constant area case is 
given at the end of the appendix. The total energy term for the stage is calculated the same as 
before (lines 20 to 31). As shown in appendix B, a total momentum source term for the stage 
can also be calculated. The actual procedure is given on lines 33 to 65. The loop for calculating 
the source terms at each grid point (lines 67 to 98) is reduced to calculating the grid cell length 
(lines 69 to 94) and distributing the total source terms to the grid point in proportion to the cell 
length (lines 95 to 97). The simplification due to constant area is evident from the much shorter 
DO 740 loop, compared to the loop for the variable area version. 
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1 C **> HEAT RELEASE section of SUBROUTINE USS as modified for 

2 C **> testing of COMPRESSION (FAN) STAGE SOURCE TERMS 

3 C **> (Version used for variable flow area through stage) 

4 C **> 

5 IF(NHEAT .LE. 0) GOTO 800 

6 DO 710 I=1,NHEAT 

7 JBEG=IHTST(I) 

* 8 JEND=IHTEND(I) 

9 XHCELF = 0.5E+0*(X(JBEG)+X(JBEG+1)) 

10 IF (XHEAT(I) .GT. XHCELF) JSTAGE=JBEG 

11 IF (XHEAT(I) .LE. XHCELF) JSTAGE=JBEG-1 

12 NPTS=JEND - JBEG + 3 

13 DO 720 LC=1,20 

14 720 HEATN(LC)=HEATIN(LC,I) 

15 TIME=TAU /TIMF AC 

16 C Get the current total pressure ratio HEATT(I) 

17 CALL LININT(20,TIME,HEATT(I),TIMHT,HEATN,IDUM) 

18 IF (HEATT(I) .LE. 1.0) GO TO 710 

19 PR=HEATT(I) 

20 ETA=XHEAT(2) 

21 C 

22 C If using a compressor map 

23 C (1) calculate corrected airflow at JSTAGE 

24 C (2) interpolate for mechanical speed (external input) and calculate 

25 C corrected speed 

26 C (3) call map subroutine to get total pressure ratio (PR) and efficiency 

27 C (ETA) from performance map 

28 C 

29 C Calculate the overall temperature delta (TT3-TT2)/TT0 

30 DTO VRT=(PR* * (GM1 / G AM)-1 .0) /ETA 

31 UXX=Q( JSTAGE, 2) / Q(JSTAGE,1) 

32 EXX=Q( JSTAGE, 3) /Q( JSTAGE, 1) 

33 C Constant in TTXX (1.2) & PTXX (1.35209) calcs assumes GAMA=1.4 

34 TTJSTAGE=GM1*(GAM*EXX - 0.5*GM1*UXX*UXX)/1.2 

35 DTOVRT=DTOVRT*TTJSTAGE 

36 CELHRT3=3.0*DTOVRT 

37 C Calculate the overall pressure delta (PT3-PT2)/PTO 

38 PXX=GM1*(Q(JSTAGE,3)-.5*Q(JSTAGE,2)*UXX) 

39 TXX=GAM*PXX/Q(JSTAGE,1) 

40 PXX=PXX*PSX(JSTAGE)/AR(JSTAGE) 

41 PTJSTAGE=PXX*((1.2*TTJSTAGE/TXX)**(GAM/GM1))/1.35209 

42 DPOVRP=(PR - 1 .0) * PT J STAGE 

43 C Begin loop for calculating source terms at each grid point 

44 DO 740 J=JBEG,JEND 

45 DXCEL=0.5E-f0*(X(J+l)-X(J-l)) 

46 IF(J .NE. JBEG) GO TO 732 
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47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 


XHCELF = 0.5E+0*(X(J)+X(J+1)) 

IF(XHEAT(I) .GT. XHCELF) DXCEL=0.0E+0 
C**> 

C**> IF DXCEL=0 the source terms will be 0 ==> go to end of DO LOOP 

e**> 

IF(DXCEL .EQ. 0.) GO TO 740 

IF(XHEAT(I) .LE. XHCELF) DXCEL=0.5E+0*(X(J)+X(J+1))-XHEAT(I) 
732 CONTINUE 

IF(J .NE. JBEG+1) GO TO 734 
XHCELB = 0.5E+0*(X(J)+X(J-1)) 

IF(XHEAT(I) .GT. XHCELB) DXCEL=0.5E44)*(X(J)+X(J+1))-XHEAT(I) 
734 CONTINUE 

IF(J .NE. JEND-1) GO TO 736 
XHCELF = 0.5E+0*(X(J)+X(J+1)) 

IF(ENDHT(I) .LT. XHCELF) DXCEL=ENDHT(I)- 
1 0.5E+0*(X(J)+X(J-1)) 

736 CONTINUE 

IF(J .NE. JEND) GO TO 738 
XHCELB = 0.5E+0*(X(J)+X(J-1)) 

IF(ENDHT(I) .LE. XHCELB) DXCEL=0.0E+0 
IF(ENDHT(I) .GT. XHCELB) DXCEL=ENDHT(I)- 
1 0.5E+0*(X(J)+X(J-1)) 

738 CONTINUE 
C 

C Calculate the cell work term for the energy equation 
C 

CELHR=(DXCEL/LHEAT(I))*CELHRT3 

G(J,3)=G(J,3)+CELHR*Q(J,2)*PSX(J) 

C 

C Calculate the cell force term for the momentum equation 


C START local avg (J) variables 

AR1=(AR(J) + AR(J-l))/2.0 
AR2=(AR(J) + AR(J+l))/2.0 
if (j .ne. jend) go to 760 

CALL LININT(NPTS,endht(i),ar2,x(jbegml),ar(jbegml),IDUM) 
760 CONTINUE 

Q1AVG=(Q(J-1,1) + Q(J,l))/2.0 
Q2AVG=(Q(J-1,2) + Q(J,2))/2.0 
Q3AVG=(Q(J-1,3) + Q(J,3))/2.0 
UXX=Q2AVG/Q1AVG 
EXX=Q3AVG/Q1AVG 

TTXX=GM1*(GAM*EXX - 0.5*GM1*UXX*UXX) /1.2 

PXX=GM1*(Q3AVG-.5*Q2AVG*UXX) 

TXX=GAM*PXX/Q1AVG 
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93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 

117 

118 

119 

120 
121 
122 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 


XXM=UXX/SQRT(TXX) 

PXX=P XX* P SX( J ) / AR 1 

PTXX=PXX*((1.2*TTXX/TXX)**(GAM/GM1))/1. 35209 
C END local avg (J) variables 


if (j .ne. jbeg) go to 759 

CALL LININT(NPTS,xheat(i) ,qlavg,x(jbegml) ,q(jbegml,l) ,IDUM) 
CALL LININT(NPTS ,xheat (i) ,q2avg,x (j begml ) ,q(j begml ,2} ,ID UM) 
CALL LININT(NPTS,xheat(i),q3avg,x(jbegml),q(jbegml,3),IDUM) 
CALL LININT(NPTS,xheat(i),arl,x(jbegml),ar(jbegml),IDUM) 
UXX=Q2AVG/Q1AVG 
EXX=Q3AVG/Q1AVG 

TTXX=GM1 * (GAM* EXX - 0.5*GM1*UXX*UXX)/1.2 
PXX=GM1 * (Q3 A VG- . 5 * Q2 A VG * UXX) 
TXX=GAM*PXX/Q1AVG 
XXM=UXX/SQRT(TXX) 

PTXX=PTJSTAGE 
TTXX=TT J S TAG E 
759 CONTINUE 


C Calculate tot. temp, and press, ratio for cell 

TTR=DTOVRT*DXCEL/TTXX/LHEAT(I) + 1.0 
PTR=DPOVRP*DXCEL/PTXX/LHEAT(I) + 1.0 
C 

C Iterate for cell exit Mach no. 

C CAUTION! If CF2>CFMAX the iteration will not converge, 
C (CFMAX=CF(M=1.0)=1.2**-3=.5787037) 

C 

CFl=XXM*(1.0 + 0.2*XXM*XXM)**-3 
CF2=ARl*CFl*SQRT(TTR)/(PTR*AR2) 

IF (XXM .GT. 1.0) GO TO 751 
C Subsonic cell inlet Mach no. 

XMGUESS=0.3 

753 XXM2=CF2*(1.0 + 0.2*XMGUESS*XMGUESS)**3 
DELM=ABS(XMGUESS-XXM2) 

IF (DELM .LT. 0.0001) GO TO 752 

XMGUESS=XXM2 

GO TO 753 

C Supersonic cell inlet Mach no. 

751 XMGUESS=1.5 

754 XXM2=5.0*((XMGUESS/CF2)**0.33333 - 1.0) 
XXM2=SQRT(XXM2) 

DELM=ABS(XMGUESS-XXM2) 

IF (DELM .LT. 0.0001) GO TO 752 
XMGUESS=XXM2 
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139 


GO TO 754 

140 

752 

CONTINUE 

141 


XMRAT=(1.0 + 0.2*XXM*XXM)/(1.0 + 0.2*XXM2*XXM2) 

142 


FAC1=SQRT(XMRAT)*XXM2*SQRT(TTR)/XXM - 1.0 

143 


FAC2=(AR1 +AR2)*(PTR*XMRAT**3.5 - 1.0)/(2.Q*ARl) 

144 


FFAC1=Q2AVG*Q2AVG*FAC1/Q1AVG 

145 


FF AC2=GM1*F AC2*(Q3AV G - 0.5*Q2AVG*Q2AVG/Q1AVG) 

146 


FFAC=(FFAC1 + FFAC2) 

147 


G(J,2) = G(J,2) + PSX(J)*FFAC 

148 

740 

CONTINUE 

149 

710 

CONTINUE 

150 

800 

CONTINUE 
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1 C **> HEAT RELEASE section of SUBROUTINE USS as modified for 

2 C **> testing of COMPRESSION (PAN) STAGE SOURCE TERMS 

3 C **> (Version used for constant flow area through stage) 

4 C **> 

5 IF(NHEAT .LE. 0) GOTO 800 

6 DO 710 I=1,NHEAT 

7 DO 720 LC=1,20 

8 720 HEATN(LC)=HEATIN(LC,I) 

9 TIME=TAU/TIMFAC 

10 CALL LININT(20,TIME,HEATT(I),TIMHT,HEATN,IDUM) 

11 JBEG=IHTST(I) 

12 JEND=IHTEND(I) 

13 

14 C **> Begin additions/changes 

15 XHCELF = 0.5E+0*(X(JBEG)+X(JBEG+1)) 

16 IF (XHEAT(I) .GT. XHCELF) J STAGE= JBEG 

17 IF (XHEAT(I) .LE. XHCELF) JSTAGE=JBEG-1 

18 IF (HEATT(I) .LE. 1.0) GO TO 710 

19 C 

20 C Calculate the work term for the energy equation 

21 C 

22 C HEATT(I) is fan PR & XHEAT(2) is used as the fan efficiency 

23 C 

24 DTOVRT=(HEATT(I)**(GM1/GAM)-1.0)/XHEAT(2) 

25 TTR=DTOVRT + 1.0 

26 UXX=Q(JSTAGE,2)/Q(JSTAGE,1) 

27 EXX=Q(JSTAGE,3)/Q(JSTAGE,1) 

28 C Constant in TTXX (1.2) calc assumes GAMA=1.4 

29 TTXX=GM1*(GAM*EXX - 0.5*GM1*UXX*UXX)/1.2 

30 DTOVRT=DTOVRT*TTXX 

31 CELHRT=3.0*DTOVRT/LHEAT(I) 

32 C 

33 C Calculate the total force term for the momentum equation 

34 C 

35 PXX=GM1*(Q(JSTAGE,3)-.5*Q(JSTAGE,2)*UXX) 

36 TXX=GAM*PXX/Q(JSTAGE,1) 

37 XXM=UXX/SQRT(TXX) 

38 CF1=XXM*(1.0 + 0.2*XXM*XXM)**-3 

39 CF2=CF1*SQRT(TTR)/HEATT(I) 

40 C**> 

41 C**> CFMAX=(1.2)**-3=.5787037 

42 C**> IF CF2>CFMAX The following iteration will not converge! 

43 C**> (not a likely situation - CF2 should be < CFl) 

44 C**> 

45 IF (XXM .GT. 1.0) GO TO 751 

46 XMGUESS=0.3 
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47 753 XXM2=CF2*(1.0 + 0.2*XMGUESS*XMGUESS)**3 

48 DELM=ABS(XMGUESS-XXM2) 

49 IF (DELM .LT. 0.0001) GO TO 752 

50 XMGUESS=XXM2 

51 GO TO 753 

52 751 XMGUESS=1.5 

53 754 XXM2=5.0*((XMGUESS/CF2)**0.33333 - 1.0) 

54 XXM2=SQRT(XXM2) 

55 DELM=ABS(XMGUESS-XXM2) 

56 IF (DELM .LT. 0.0001) GO TO 752 

57 XMGUESS=XXM2 

58 GO TO 754 

59 752 CONTINUE 

60 XMRAT=(1.0 + 0.2*XXM*XXM)/(1.0 + 0.2*XXM2*XXM2) 

61 FAC1=SQRT(XMRAT)*XXM2*SQRT(TTR)/XXM - 1.0 

62 FAC2=HEATT(I)*XMRAT**3.5 - 1.0 

63 FFAC1=Q(JSTAGE,2)*Q(JSTAGE,2)*FAC1/Q(JSTAGE,1) 

64 FFAC2=PXX*FAC2 

65 FFAC=FFAC1 + FFAC2 

66 C **> End additions/changes 

67 C Begin loop for calculating source terms at each grid point 

68 DO 740 J=JBEG,JEND 

69 DXCEL=0.5E+0*(X(J+1)-X(J-1)) 

70 IF(J .NE. JBEG) GO TO 732 

71 XHCELF = 0.5E+0*(X(J)+X(J+1)) 

72 C**> 

73 IF(XHEAT(I) .GT. XHCELF) DXCEL=0.0E+0 

74 C**> 

75 C**> If DXCEL=0 the source terms are 0 ==> go to end of DO LOOP 

76 C**> 

77 IF(DXCEL .EQ. 0.) GO TO 740 

78 IF(XHEAT(I) .LE. XHCELF) DXCEL=0.5E+0*(X(J)+X(J+l))-XHEAT(I) 

79 732 CONTINUE 

80 IF(J .NE. JBEG+1) GO TO 734 

81 XHCELB = 0.5E+0* (X( J) +X( J- 1)) 

82 IF(XHEAT(I) .GT. XHCELB) DXCEL=0.5E+0*(X(J)+X(J+1))-XHEAT(I) 

83 734 CONTINUE 

84 IF(J .NE. JEND-1) GO TO 736 

85 XHCELF = 0.5E+0*(X(J)+X(J+1)) 

86 IF(ENDHT(I) .LT. XHCELF) DXCEL— ENDHT(I)- 

87 1 Q.5E+0*(X(J)+X(J-1)) 

88 736 CONTINUE 

89 IF(J .NE. JEND) GO TO 738 

90 XHCELB = 0.5E+0*(X(J)+X(J-1)) 

91 IF(ENDHT(I) .LE. XHCELB) DXCEL=0.0E+0 

92 IF(ENDHT(I) .GT. XHCELB) DXCEL=ENDHT(I)- 
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0.5E+0*(X(J)+X(J-1)) 


93 

] 

l 

94 

738 

CONTINUE 

95 


CELHR=DX< 

96 


C(J,3)=G(J,3 

97 


G(J,2) = G(J 

98 

740 

CONTINUE 

99 

710 

CONTINUE 

100 

800 

CONTINUE 
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Figure 3. — Compression stage characteristics. 
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***** INPUT DATA MIRROR ***** 


***************** SIMPLE GEOMETRY WITH NO CENTERBODY **************** 

************************ Data for demonstration case ************************ 
SCONTRL SCMC=.TRUE., NMAX=2800, NC=20, TIMMAX=.30, 

TlMDTC=0.0D+0, .025, DT1M= .0001, .0001, NDTC=2 

SEND 

SBC HHUSIC=2.0, MHDSBC=.TRUE., MHDSIC=.164, 

PRUSIC= 270.42, TMUSIC=288.36, 

TIMDS=0.0D+0, .09, .12, 100.0, MACHDS= 2*1.0, 2*2.0622, 

START=.TRUE., OPCOND*. FALSE., 

SEND 

SGEOM XC=0. 00+0,. 7, 1.4, 2. 5, 3. 5, 

YC=1 .0, .81 ,1.0,1.05,1.5, 

XS=0.,3.6, 

YS=2*0. 

TWOO=. FALSE., INTYPE*0, XUSLIM=.0, 

XINLET=0., XEXIT=3.5, JMAX=36, JINL=1, 

IPTRAN=8, 

COULR- .8333333 

SEND 

SBYBL NBL=0 SEND 
SINJECT SEND 
SHEAT NHEAT-1, 

XHEAT=1.56, .65, LHEAT=.7, 

TIMHT=0.0D+0, 0.03, 0.07, .14, .24, 100.0, 

HEATIN* 2*1.0, 2*1.4, 2*2.4, 

SEND 

SLOSS SEND 

***** END DATA INPUT ***** 

***************** SIMPLE GEOMETRY WITH NO CENTERBOOY **************** 

************************ Data for demonstration case ************************ 


IMPLICIT SPLIT CHARACTERISTICS 

AXISYMMETRIC - VARIABLES NORMALIZED BY INITIAL UPSTREAM STATIC OR TOTAL CONDTIONS 
UPSTREAM MACH NUMBER BOUNDARY CONDITION, MHUSIC = 2.000000 

DOWNSTREAM MACH NO. BOUNDARY CONDITION, MHDSIC = 0.164000 

N= 300 TIME* 0.0300 SEC XSHK* 0.0000 XSH2* 0.0000 


X 

RHO 

U 

T 

P 

M 

CEL 

MASS 

AREA 

FT 

TT 

CMASS 

0.0000 

1.0000 

1.0000 

1.0000 

1.0000 

2.0000 

2.9968 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

0.1000 

1.0771 

0.9810 

1.0301 

1.1095 

1.9330 

2.9738 

1.0000 

0.9465 

1.0000 

1.0000 

1.0000 

0.2000 

1.1649 

0.9598 

1.0630 

1.2382 

1.8620 

2.9476 

1.0000 

0.8944 

1.0000 

1.0000 

1.0000 

0.3000 

1.2661 

0.9361 

1.0990 

1.3914 

1.7859 

2.9174 

1.0000 

0.8438 

1.0000 

1.0000 

1.0000 

0,4000 

1.3844 

0.9090 

1.1390 

1.5768 

1.7034 

2.8822 

1.0000 

0.7946 

0.9999 

1.0000 

1.0001 

0.5000 

1.5259 

0.8774 

1.1842 

1.8069 

1.6125 

2.8399 

1.0000 

0.7470 

0.9999 

1.0000 

1.0001 

0.6000 

1.7005 

0.8391 

1.2367 

2.1031 

1.5091 

2.7873 

1.0000 

0.7008 

0.9998 

1.0000 

1.0002 

0.7000 

1.7773 

0.8404 

1.2485 

2.2190 

1.5042 

2.7951 

0.9799 

0.6561 

1.0474 

1.0075 

0.9391 

0.8000 

3.2229 

0.4447 

1.6422 

5.2927 

0.6941 

2.1687 

1.0045 

0.7008 

0.9334 

1.0003 

1.0764 

0.9000 

3.4609 

0.3875 

1.6764 

5.8019 

0.5985 

2.0675 

1.0017 

0.7470 

0.9447 

0.9981 

1.0593 

1.0000 

3.5584 

0.3536 

1.7002 

6.0500 

0.5424 

2.0090 

0.9999 

0.7946 

0.9445 

1.0001 

1.0588 

1,1000 

3.6403 

0.3258 

1.7146 

6.2415 

0.4976 

1.9589 

1.0007 

0.8438 

0.9447 

0.9997 

1.0591 

1.2000 

3.7017 

0.3021 

1.7274 

6.3944 

0.4597 

1.9165 

1.0002 

0.8944 

0.9446 

1.0002 

1.0589 

1.3000 

3.7542 

0.2815 

1.7365 

6.5193 

0.4273 

1.8789 

1.0004 

0.9465 

0.9446 

1.0000 

1 .0590 

1.4000 

3.7964 

0.2635 

1.7446 

6.6231 

0.3989 

1.8458 

1.0002 

1.0000 

0.9446 

1,0001 

1.0589 

1.5000 

3.8028 

0.2606 

1.7457 

6.6386 

0.3945 

1.8405 

1.0001 

1.0091 

0.9445 

1.0000 

1.0589 

1 .6000 

3.8089 

0.2578 

1.7469 

6.6538 

0.3902 

1.8354 

1.0000 

1.0183 

0.9445 

1.0000 

1.0588 

1.7000 

3.8150 

0.2551 

1.7480 

6.6685 

0.3859 

1.8304 

1.0000 

1.0275 

0.9445 

1.0000 

1.0588 

1.8000 

3.8209 

0.2525 

1.7490 

6.6829 

0.3818 

1.8255 

1.0000 

1.0367 

0.9445 

1.0000 

1.0588 

1.9000 

3.8267 

0.2498 

1.7500 

6.6969 

0.3777 

1.8207 

1.0000 

1.0460 

0.9445 

1.0000 

1.0588 

2.0000 

3.8323 

0.2473 

1.7511 

6.7106 

0.3737 

1.8159 

1.0000 

1.0553 

0.9445 

1.0000 

1.0588 

2.1000 

3.8378 

0.2448 

1.7520 

6.7239 

0.3698 

1.8112 

1.0000 

1.0646 

0.9445 

1.0000 

1.0588 

2.2000 

3.8430 

0.2423 

1.7530 

6.7369 

0.3660 

1.8067 

1.0000 

1 .0740 

0.9445 

1.0000 

1.0588 

2.3000 

3.8482 

0.2398 

1.7540 

6.7495 

0.3622 

1.8022 

1.0000 

1.0835 

0.9445 

1.0000 

1.0588 

2.4000 

3.8531 

0.2375 

1.7549 

6.7619 

0.3585 

1.7977 

1.0000 

1.0930 

0.9445 

1.0000 

1.0588 

2.5000 

3.8580 

0.2351 

1.7558 

6.7738 

0.3549 

1.7934 

1.0000 

1.1025 

0.9445 

1.0000 

1.0588 

2.6000 

3.9001 

0.2138 

1.7634 

6.8775 

0.3221 

1.7538 

1.0000 

1.1990 

0.9445 

1.0000 

1.0588 

2.7000 

3.9331 

0.1956 

1.7694 

6.9592 

0.2941 

1.7196 

1.0000 

1.2996 

0.9445 

1.0000 

1.0588 

2.8000 

3.9595 

0.1798 

1.7741 

7.0247 

0.2700 

1.6899 

1.0000 

1.4042 

0.9445 

1.0000 

1.0588 

2.9000 

3.9809 

0.1660 

1.7779 

7.0778 

0.2490 

1.6637 

0.9999 

1.5129 

0.9445 

1.0000 

1.0588 

3.0000 

3.9984 

0.1538 

1.7811 

7.1213 

0.2305 

1.6405 

0.9999 

1.6256 

0.9445 

1.0000 

1.0588 

3.1000 

4.0128 

0.1430 

1.7836 

7.1574 

0.2142 

1.6198 

0.9999 

1.7424 

0.9445 

1.0000 

1.0588 


Figure 4. — Listing of input data and results at selected time steps for simulated fan starting transient. 
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3.2000 

4.0249 

0.1333 

1.7858 

7.1874 

0.1996 

1.6013 

1.0000 

1.8632 

0.9445 

1.0000 

1.0588 

3.3000 

4.0350 

0.1247 

1.7876 

7.2127 

0.1865 

1.5846 

1.0000 

1.9881 

0.9445 

1.0000 

1.0588 

3.4000 

4.0435 

0.1168 

1.7891 

7.2341 

0.1747 

1.5695 

1.0000 

2.1170 

0.9445 

1.0000 

1.0588 

3.5000 
1000 TIME 

4.0508 0.1097 1.7903 7.2524 

[= 0.1000 SEC XSHK = 0.0000 XSH2 = 

0.1640 

0.0000 

1.5558 

1.0000 

2.2500 

0.9445 

1.0000 

1.0588 

X 

RHO 

U 

T 

P 

H 

CFl 

MASS 

AREA 

PT 

TT 

CMASS 

0.0000 

1.0000 

1.0000 

1.0000 

1.0000 

2.0000 

2.9968 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

0.1000 

1.0771 

0.9810 

1.0301 

1.1095 

1.9330 

2.9738 

1.0000 

0.9465 

1.0000 

1.0000 

1.0000 

0.2000 

1.1649 

0.9598 

1.0630 

1.2382 

1.8620 

2.9476 

1.0000 

0.8944 

1.0000 

1.0000 

1.0000 

0.3000 

1.2661 

0.9361 

1.0990 

1.3914 

1.7859 

2.9174 

1.0000 

0.8438 

1.0000 

1.0000 

1.0000 

0.4000 

1.3844 

0.9090 

1.1390 

1.5768 

1.7034 

2.8822 

1.0000 

0.7946 

0.9999 

1.0000 

1.0001 

0.5000 

1.5259 

0.8774 

1.1842 

1.8069 

1.6125 

2.8399 

1.0000 

0.7470 

0.9999 

1.0000 

1.0001 

0.6000 

1.7005 

0.8391 

1.2367 

2.1031 

1.5091 

2.7873 

1.0000 

0.7008 

0.9998 

1.0000 

1.0002 

0.7000 

1.9299 

0.7898 

1.3010 

2.5109 

1.3848 

2.7172 

1.0000 

0.6561 

0.9996 

1.0000 

1.0004 

0.8000 

1.7005 

0.8391 

1.2367 

2.1031 

1.5091 

2.7873 

1.0000 

0.7008 

0.9998 

1.0000 

1.0002 

0.9000 

1.5259 

0.8774 

1.1842 

1.8069 

1.6125 

2.8399 

1.0000 

0.7470 

0.9999 

1.0000 

1.0001 

1.0000 

1.3844 

0.9090 

1.1390 

1.5768 

1.7034 

2.8822 

1.0000 

0.7946 

0.9999 

1.0000 

1.0001 

1.1000 

1.2661 

0.9361 

1.0990 

1.3914 

1.7859 

2.9174 

1.0000 

0.8438 

1.0000 

1.0000 

1.0000 

1.2000 

1.1649 

0.9598 

1.0630 

1.2382 

1.8620 

2.9476 

1.0000 

0.8944 

1.0000 

1.0000 

1.0000 

1.3000 

1.0771 

0.9810 

1.0301 

1.1095 

1.9330 

2.9738 

1.0000 

0.9465 

1.0000 

1.0000 

1.0000 

1.4000 

1.0000 

1.0000 

1.0000 

1.0000 

2.0000 

2.9968 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.5000 

0.9880 

1.0030 

0.9952 

0.9832 

2.0109 

3.0004 

1.0000 

1.0091 

1.0000 

1.0000 

1.0000 

1.6000 

0.9654 

1.0173 

0.9901 

0.9558 

2.0448 

3.0264 

1.0000 

1.0183 

1.0248 

1.0100 

0.9807 

1.7000 

0.9327 

1.0435 

0.9848 

0.9186 

2.1030 

3.0761 

1.0000 

1.0275 

1.0785 

1.0310 

0.9415 

1.8000 

0.9018 

1.0696 

0.9805 

0.8843 

2.1604 

3.1261 

1.0000 

1.0367 

1.1358 

1.0532 

0.9036 

1.9000 

0.8735 

1.0945 

0.9773 

0.8537 

2.2143 

3.1743 

1.0000 

1.0460 

1.1931 

1.0754 

0.8692 

2.0000 

1 .9769 

0.5512 

1.6686 

3.2986 

0.8535 

2.3917 

1.1499 

1.0553 

0.6786 

1.0620 

1.7464 

2.1000 

2.7080 

0.3561 

1.9046 

5.1576 

0.5160 

2.0900 

1.0265 

1.0646 

0.7904 

1.1145 

1.3710 

2.2000 

2.8922 

0.3317 

1.9573 

5.6611 

0.4741 

2.0602 

1.0303 

1.0740 

0.8439 

1.1363 

1.3014 

2.3000 

2.9890 

0.3194 

1.9844 

5.9312 

0.4535 

2.0454 

1.0345 

1.0835 

0.8729 

1.1478 

1.2697 

2.4000 

2.9956 

0.3174 

1.9851 

5.9464 

0.4506 

2.0417 

1.0393 

1.0930 

0.8736 

1.1476 

1.2745 

2.5000 

2.9958 

0.3165 

1.9814 

5.9359 

0.4497 

2.0385 

1.0453 

1.1025 

0.8716 

1.1453 

1.2836 

2.6000 

3.0518 

0.2878 

1.9883 

6.0681 

0.4082 

1.9836 

1.0532 

1.1990 

0.8698 

1.1414 

1.2936 

2.7000 

3.1030 

0.2637 

1.9883 

6.1695 

0.3740 

1.9354 

1.0633 

1.2996 

0.8684 

1.1355 

1.3047 

2.8000 

3.1531 

0.2430 

1.9819 

6.2491 

0.3453 

1.8919 

1.0761 

1.4042 

0.8673 

1.1273 

1.3173 

2.9000 

3.2037 

0.2251 

1.9705 

6.3129 

0.3208 

1.8520 

1.0912 

1.5129 

0.8664 

1.1172 

1.3312 

3.0000 

3.2540 

0.2095 

1.9560 

6.3649 

0.2995 

1.8156 

1.1080 

1.6256 

0.8657 

1.1062 

1.3461 

3.1000 

3.3023 

0.1956 

1.9405 

6.4080 

0.2808 

1.7823 

1.1254 

1.7424 

0.8651 

1.0950 

1.3614 

3.2000 

3.3460 

0.1833 

1.9261 

6.4447 

0.2641 

1.7525 

1.1425 

1.8632 

0.8646 

1.0850 

1.3764 

3.3000 

3.3820 

0.1722 

1.9149 

6.4762 

0.2489 

1.7264 

1.1579 

1.9881 

0.8641 

1.0770 

1.3906 

3.4000 

3.4073 

0.1623 

1.9087 

6.5034 

0.2349 

1 .7043 

1.1707 

2.1170 

0.8637 

1.0721 

1.4034 

3.5000 
2800 TIME ; 

3.4202 0.1533 1.9085 6.5273 

= 0.2800 SEC XSHK = 0.0000 XSH2 = 

0.2220 

0.0000 

1.6864 

1.1801 

2.2500 

0.8634 

1.0707 

1.4143 

X 

RHO 

U 

T 

P 

M 

CFl 

MASS 

AREA 

PT 

TT 

CMASS 

0.0000 

1.0000 

1.0000 

1.0000 

1.0000 

2.0000 

2.9968 

1.0000 

1.0000 

1.0000 

1.0000 

1 .0000 

0.1000 

1.0771 

0.9810 

1.0301 

1.1095 

1.9330 

2.9738 

1.0000 

0.9465 

1.0000 

1.0000 

1 .0000 

0.2000 

1.1649 

0.9598 

1.0630 

1.2382 

1.8620 

2.9476 

1.0000 

0.8944 

1.0000 

1.0000 

1.0000 

0.3000 

1.2661 

0.9361 

1.0990 

1.3914 

1.7859 

2.9174 

1.0000 

0.8438 

1.0000 

1.0000 

1.0000 

0.4000 

1 .3844 

0.9090 

1.1390 

1.5768 

1.7034 

2.8822 

1.0000 

0.7946 

0.9999 

1.0000 

1.0001 

0.5000 

1.5259 

0.8774 

1.1842 

1.8069 

1.6125 

2.8399 

1.0000 

0.7470 

0.9999 

1.0000 

1 .0001 

0.6000 

1.7005 

0.8391 

1.2367 

2.1031 

1.5091 

2.7873 

1.0000 

0.7008 

0.9998 

1.0000 

1 .0002 

0.7000 

1.9299 

0.7898 

1.3010 

2.5109 

1.3848 

2.7172 

1.0000 

0.6561 

0.9996 

1.0000 

1.0004 

0.8000 

1.7005 

0.8391 

1.2367 

2.1031 

1 .5091 

2.7873 

1.0000 

0.7008 

0.9998 

1.0000 

1.0002 

0.9000 

1.5259 

0.8774 

1.1842 

1.8069 

1.6125 

2.8399 

1.0000 

0.7470 

0.9999 

1.0000 

1.0001 

1.0000 

1.3844 

0.9090 

1.1390 

1.5768 

1.7034 

2.8822 

1.0000 

0.7946 

0.9999 

1.0000 

1.0001 

1.1000 

1.2661 

0.9361 

1.0990 

1.3914 

1.7859 

2.9174 

1.0000 

0.8438 

1.0000 

1.0000 

1.0000 

1.2000 

1.1649 

0.9598 

1.0630 

1.2382 

1.8620 

2.9476 

1.0000 

0.8944 

1.0000 

1.0000 

1.0000 

1.3000 

1.0771 

0.9810 

1.0301 

1.1095 

1.9330 

2.9738 

1.0000 

0.9465 

1.0000 

1 .0000 

1.0000 

1.4000 

1.0000 

1.0000 

1.0000 

1 .0000 

2.0000 

2.9968 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.5000 

0.9880 

1.0030 

0.9952 

0.9832 

2.0109 

3.0004 

1.0000 

1.0091 

1.0000 

1.0000 

1.0000 

1 .6000 

0.9442 

1.0402 

0.9851 

0.9300 

2.0960 

3.0696 

1.0000 

1.0183 

1.0803 

1.0281 

0.9386 

1.7000 

0.8764 

1.1105 

0.9709 

0.8509 

2.2540 

3.2029 

1.0000 

1.0275 

1.2655 

1.0874 

0.8240 

1.8000 

0.8202 

1.1761 

0.9632 

0.7900 

2.3967 

3.3301 

1.0000 

1.0367 

1.4685 

1.1499 

0.7302 

1.9000 

0.7742 

1.2348 

0.9625 

0.7452 

2.5173 

3.4470 

1.0000 

1.0460 

1.6717 

1.2124 

0.6587 

2.0000 

0.7354 

1.2886 

0.9663 

0.7106 

2.6217 

3.5564 

1.0000 

1.0553 

1.8740 

1.2748 

0.6025 

2.1000 

0.7016 

1 .3387 

0.9734 

0.6829 

2.7138 

3.6602 

1.0000 

1.0646 

2.0758 

1.3373 

0.5571 

2.2000 

0.6717 

1.3861 

0.9826 

0.6601 

2.7965 

3.7594 

1.0000 

1.0740 

2.2772 

1.3997 

0.5195 

2.3000 

0.6535 

1.4122 

0.9859 

0.6443 

2.8445 

3.8133 

1.0000 

1.0835 

2.3917 

1.4341 

0.5007 

2.4000 

0.6461 

1.4161 

0.9827 

0.6349 

2.8571 

3.8195 

1.0000 

1 .0930 

2.4021 

1.4372 

0.4991 

2.5000 

0.6397 

1.4178 

0.9788 

0.6262 

2.8662 

3.8210 

1.0000 

1.1025 

2.4021 

1.4372 

0.4991 

2.6000 

0.5817 

1 .4339 

0.9423 

0.5481 

2.9543 

3.8343 

1.0000 

1.1990 

2.4020 

1.4372 

0.4991 

2.7000 

0.5313 

1.4484 

0.9087 

0.4828 

3.0388 

3.8460 

1.0000 

1 .2996 

2.4020 

1.4372 

0.4991 

2.8000 

0.4872 

1.4617 

0.8778 

0.4277 

3.1202 

3.8562 

1.0000 

1.4042 

2.4019 

1.4372 

0.4991 


Figure 4 .— Continued . 
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Duct area ratio. 
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TOTAL CPU TIME: 0.1120E+03 SECONDS 

Figure 4.— Concluded. 
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Normalized distance from cowl lip 

Figure 5.— Duct area distribution for simulated inlet/fan 
system. 



(a) Diffuser exit Mach number. 



Time, sec 

(b) Fan total pressure ratio. 

Figure 6.— Fan stage and diffuser exit boundary condition 
inputs for simulated fan stage starting transient. 
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(d) Total temperature. 

Figure 7.— Axial distributions at various time steps during 
simulated fan starting transient. 




Time, sec 

(b) Total temperature ratio. 

Figure 8. — Calculated fan-ratio time histories during simu- 
lated fan starting transient. 
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